In [1]:
%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import sqrt, sin, pi, var, integrate, Symbol, S, Integral
var("x y z")
n = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
Vh = -((x-1)**4 + (y-1)**4 + (z-1)**4) + 2*(x**2+y**2+z**2)
def laplace(f):
return f.diff(x, 2) + f.diff(y, 2) + f.diff(z, 2)
print "n ="
n
Out[1]:
In [6]:
print "Vh ="
Vh
Out[6]:
Let's check charge neutrality:
In [7]:
integrate(n, (x, 0, 2), (y, 0, 2), (z, 0, 2))
Out[7]:
$V_H$ is produced by the charge density $n$:
In [8]:
laplace(Vh) / (-4*pi)
Out[8]:
In [9]:
_.simplify()
Out[9]:
The Hartree energy $E_H$ is:
In [10]:
integrate(Vh*n, (x, 0, 2), (y, 0, 2), (z, 0, 2))/2
Out[10]:
We are solving a Poisson equation: $$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ with $$n(x, y, z) = {3\over\pi} ((x-1)^2 + (y-1)^2 + (z-1)^2 - 1)$$ on the domain $[0, 2] \times [0, 2] \times [0, 2]$ with periodic boundary condition. This charge density is net neutral. The solution is given by: $$ V_H({\bf x}) = \int {n({\bf y})\over |{\bf x}-{\bf y}|} d^3 y $$ where ${\bf x}=(x, y, z)$. In reciprocal space $$ V_H({\bf G}) = 4\pi {n({\bf G})\over G^2} $$ where ${\bf G}$ are the reciprocal space vectors. I am interested in the Hartree energy: $$ E_H = {1\over 2} \int {n({\bf x}) n({\bf y})\over |{\bf x}-{\bf y}|} d^3 x d^3 y = {1\over 2} \int V_H({\bf x}) n({\bf x}) d^3 x $$ In reciprocal space this becomes (after discretization): $$ E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2} $$ The ${\bf G}=0$ term is omitted, which effectively makes the charge density net neutral (and since it is already neutral, then everything is consistent).
For the test problem above, this can be evaluated analytically and one gets: $$ E_H = {128\over 35\pi} = 1.16410... $$ Now we implement the FFT solver in NumPy below and plot the convergence graphs with respect to this exact energy.
This problem has also been discussed at:
http://scicomp.stackexchange.com/questions/7097/convergence-rate-of-fft-poisson-solver/7169
Computational Science StackExchange is a great site to ask these questions!
In [11]:
from numpy import empty, pi, meshgrid, linspace, sum
from numpy.fft import fftn, fftfreq
E_exact = 128/(35*pi)
print "Hartree Energy (exact): %.15f" % E_exact
f = open("conv.txt", "w")
for N in range(3, 200, 10):
print "N =", N
L = 2.
x1d = linspace(0, L, N+1)[:-1]
x, y, z = meshgrid(x1d, x1d, x1d)
nr = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
ng = fftn(nr) / N**3
G1d = N * fftfreq(N) * 2*pi/L
kx, ky, kz = meshgrid(G1d, G1d, G1d)
G2 = kx**2+ky**2+kz**2
G2[0, 0, 0] = 1 # omit the G=0 term
tmp = 2*pi*abs(ng)**2 / G2
tmp[0, 0, 0] = 0 # omit the G=0 term
E = sum(tmp) * L**3
print "Hartree Energy (calculated): %.15f" % E
f.write("%d %.15f\n" % (N, E))
f.close()
In [12]:
from pylab import plot, semilogy
from numpy import loadtxt, pi
D = loadtxt("conv.txt")
N = D[:, 0]
E = D[:, 1]
E_exact = 128/(35*pi)
In [13]:
figure(figsize=(8, 6), dpi=80)
semilogy(N, E-E_exact, "k-")
grid()
title("Convergence of an FFT Poisson solver (semilogy)")
xlabel("N (number of PW in each direction)")
ylabel("E - E_exact [a.u.]")
savefig("fft_convergence_semilogy.png")
In [14]:
figure(figsize=(8, 6), dpi=80)
calculated = E-E_exact
loglog(N, E-E_exact, "ko", label="calculated")
predicted = 1/N**2
predicted = predicted * (calculated[-1] / predicted[-1])
loglog(N, predicted, "g-", label="$1/N^2$")
grid()
title("Convergence of an FFT Poisson solver (loglog)")
xlabel("N (number of PW in each direction)")
ylabel("E - E_exact [a.u.]")
legend()
savefig("fft_convergence_loglog.png")
Rerun the above FFT calculation with the following potential:
In [3]:
nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4
Plot convergence study. Why is the convergence rate different?
In [9]:
from numpy import empty, pi, meshgrid, linspace, sum, exp, sin
from numpy.fft import fftn, fftfreq
f = open("conv2.txt", "w")
for N in range(3, 30, 2):
print "N =", N
L = 2.
x1d = linspace(0, L, N+1)[:-1]
x, y, z = meshgrid(x1d, x1d, x1d)
nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4
ng = fftn(nr) / N**3
G1d = N * fftfreq(N) * 2*pi/L
kx, ky, kz = meshgrid(G1d, G1d, G1d)
G2 = kx**2+ky**2+kz**2
G2[0, 0, 0] = 1 # omit the G=0 term
tmp = 2*pi*abs(ng)**2 / G2
tmp[0, 0, 0] = 0 # omit the G=0 term
E = sum(tmp) * L**3
print "Hartree Energy (calculated): %.15f" % E
f.write("%d %.15f\n" % (N, E))
f.close()
In [11]:
from pylab import plot, semilogy
from numpy import loadtxt, pi
D = loadtxt("conv2.txt")
N = D[:, 0]
E = D[:, 1]
E_exact = E[-1]
figure(figsize=(8, 6), dpi=80)
calculated = E-E_exact
semilogy(N, E-E_exact, "ko-", label="calculated")
grid()
title("Convergence of an FFT Poisson solver (loglog)")
xlabel("N (number of PW in each direction)")
ylabel("E - E_exact [a.u.]")
legend()
savefig("fft_convergence_loglog.png")
In [ ]: